Check the response code. This seems to generate errors (sometimes).
# test
url <- "https://api.octopus.energy/v1/products"
message("Getting: ", url)
## Getting: https://api.octopus.energy/v1/products
resp <- httr::GET(url)
message("Status code: ", resp$status_code)
## Status code: 200
# df <- jsonlite::parse_json(resp, simplifyVector = TRUE)
# makeFlexTable(head(df$results), cap = "Example products list (first 6 rows)")
url <- paste0("https://api.octopus.energy/v1/accounts/", apiParams$accountNo , "/")
resp <- httr::GET(url = url, authenticate(user = apiParams$key, password = ""))
df <- jsonlite::parse_json(resp, simplifyVector = TRUE)
## No encoding supplied: defaulting to UTF-8.
props <- data.table::as.data.table(df$properties)
makeFlexTable(head(props[, .(town, county,
electricity_meter_points,gas_meter_points)]), cap = "Properties linked to this account (non-disclosive data)")
town | county | electricity_meter_points | gas_meter_points |
BRADFORD-ON-AVON | WILTSHIRE | [[data.frame]] | [[data.frame]] |
FRAMLINGHAM | [[data.frame]] | [[data.frame]] |
List the electricity meters by MPAN. There should be only one…
# this is a list of n mpans
length(props$electricity_meter_points)
## [1] 2
message("n MPANS listed: ", length(df$properties$electricity_meter_points))
## n MPANS listed: 2
for(n in 1:length(df$properties$electricity_meter_points)){
print(props$electricity_meter_points[n])
}
## [[1]]
## mpan profile_class consumption_standard
## 1 2000006198482 1 2853
## meters
## 1 L78C64517, 01, STANDARD, TRUE
## agreements
## 1 E-1R-SUPER-GREEN-12M-20-09-22-H, E-1R-SUPER-GREEN-12M-20-09-22-H, 2020-11-01T00:00:00Z, 2020-11-21T00:00:00Z, 2020-11-01T00:00:00Z, 2021-06-30T00:00:00+01:00
## is_export
## 1 FALSE
##
## [[1]]
## mpan profile_class consumption_standard
## 1 1050001805886 1 3675
## meters
## 1 19L3027004, 1, STANDARD, TRUE
## agreements
## 1 E-1R-SUPER-GREEN-12M-20-09-22-A, E-1R-LOYAL-FIX-12M-21-10-07-A, 2021-07-01T00:00:00+01:00, 2021-11-21T00:00:00Z, 2021-11-21T00:00:00Z, 2022-11-21T00:00:00Z
## is_export
## 1 FALSE
List the gas meters by MPRN. There should be only one…
length(df$properties$gas_meter_points)
## [1] 2
message("n MPRNS listed: ", length(df$properties$gas_meter_points))
## n MPRNS listed: 2
for(n in 1:length(df$properties$gas_meter_points)){
print(df$properties$gas_meter_points[n])
}
## [[1]]
## mprn consumption_standard meters
## 1 4256845702 18791 G4A01559730801
## agreements
## 1 G-1R-SUPER-GREEN-12M-20-09-22-H, G-1R-SUPER-GREEN-12M-20-09-22-H, 2020-11-01T00:00:00Z, 2020-11-21T00:00:00Z, 2020-11-01T00:00:00Z, 2021-06-30T00:00:00+01:00
##
## [[1]]
## mprn consumption_standard meters
## 1 7825700304 15129 E6S12725512161, E6S17944211961, NOTINSTALLED
## agreements
## 1 G-1R-LOYAL-FIX-12M-21-10-07-A, 2021-10-04T00:00:00+01:00, 2022-10-04T00:00:00+01:00
See: https://www.guylipman.com/octopus/api_guide.html#s3
Start data extraction from 1st Jan 2022 as we had smart meter (re)installed in January.
url <- paste0("https://api.octopus.energy/v1/electricity-meter-points/",
apiParams$elec_import_mpan , "/",
"meters/",
apiParams$elec_import_serial, "/",
"consumption/",
"?period_from=2022-01-01T00:00Z",
"&page_size=10000")
# get data via httr ----
resp <- httr::GET(url = url, authenticate(user = apiParams$key, password = ""))
df <- jsonlite::parse_json(resp, simplifyVector = TRUE) # creates a df of which 'results' = the data
## No encoding supplied: defaulting to UTF-8.
elecCons_dt <- data.table::as.data.table(df$results) # convert to dt
# derived variables ----
elecCons_dt <- makeDerivedVars(elecCons_dt)
maxTime <- max(elecCons_dt$dv_start)
hoursAgo <- now() - maxTime
# meter is SMETS2
elecCons_dt[, consumption_kWh := consumption] # for clarity - see https://developer.octopus.energy/docs/api/#list-consumption-for-a-meter
The data used here is up to 2022-09-19 23:30:00, which is 10.1 hours ago. In general the Octopus API seems to have data up to midnight last night.
t <- elecCons_dt[, .(nDays = uniqueN(dv_date),
sumkWh = sum(consumption_kWh),
hh_meankWh = mean(consumption_kWh),
hh_minkWh = min(consumption_kWh),
hh_maxkWh = max(consumption_kWh)),
keyby = .(month = lubridate::month(dv_date))]
makeFlexTable(t, cap = "Monthly stats", digits = 2)
month | nDays | sumkWh | hh_meankWh | hh_minkWh | hh_maxkWh |
2.00 | 6 | 65.71 | 0.26 | 0.00 | 1.32 |
3.00 | 31 | 321.52 | 0.22 | 0.00 | 1.79 |
4.00 | 30 | 275.87 | 0.19 | 0.00 | 1.43 |
5.00 | 31 | 250.11 | 0.17 | 0.00 | 1.58 |
6.00 | 30 | 251.09 | 0.17 | 0.00 | 1.25 |
7.00 | 31 | 300.75 | 0.20 | 0.00 | 1.39 |
8.00 | 31 | 312.98 | 0.21 | 0.00 | 1.52 |
9.00 | 19 | 173.23 | 0.19 | 0.00 | 1.57 |
message("Total elec to date")
## Total elec to date
sum(t$sumkWh)
## [1] 1951.277
t <- elecCons_dt[, .(nDays = uniqueN(dv_date),
sumkWh = sum(consumption_kWh)),
keyby = .(dv_date)]
message("########")
## ########
elec_kWh <- 0.24
message("Projected elec annual total @ ", elec_kWh,"p")
## Projected elec annual total @ 0.24p
annualTotal <- mean(t$sumkWh)*365
annualTotal
## [1] 3407.733
message("Projected annual elec cost @ ", elec_kWh,"p - including standing charge")
## Projected annual elec cost @ 0.24p - including standing charge
annualCost <- (annualTotal*elec_kWh)+(365*0.2401) # standing charge
annualCost
## [1] 905.4923
message("Projected mean monthly elec cost @ ", elec_kWh,"p - including standing charge")
## Projected mean monthly elec cost @ 0.24p - including standing charge
annualCost/12
## [1] 75.45769
message("########")
## ########
elec_kWh <- 0.34
message("Projected elec annual total @ ", elec_kWh,"p")
## Projected elec annual total @ 0.34p
annualTotal <- mean(t$sumkWh)*365
annualTotal
## [1] 3407.733
message("Projected annual elec cost @ ", elec_kWh,"p - including standing charge")
## Projected annual elec cost @ 0.34p - including standing charge
annualCost <- (annualTotal*elec_kWh)+(365*0.2401) # standing charge
annualCost
## [1] 1246.266
message("Projected mean monthly elec cost @ ", elec_kWh,"p - including standing charge")
## Projected mean monthly elec cost @ 0.34p - including standing charge
annualCost/12
## [1] 103.8555
TO DO: fix £ analysis to use tariff from API
URL will be something like
Figure 4.1 shows half-hourly electricity import (‘consumption’) for the current year. Spot the power cuts…
To do: mark weekends somehow
ggplot2::ggplot(elecCons_dt, aes(x = dv_date, y = dv_hms, fill = consumption_kWh)) +
geom_tile() +
theme(legend.position = "bottom") +
scale_fill_viridis_c(name = "Electricity import (kWh)") +
labs(x = "Date",
y = "Half-hour")
Figure 4.1: Half hourly electricity import (current year)
Repeat but with just the last 7 days of data - useful for checking recent appliance use and offspring effects.
today <- lubridate::today()
plotDT <- elecCons_dt[dv_date >= today - 7]
p <- ggplot2::ggplot(plotDT, aes(x = dv_date, y = dv_hms, fill = consumption_kWh)) +
geom_tile() +
theme(legend.position = "bottom") +
scale_fill_viridis_c(name = "Electricity import (kWh)") +
labs(x = "Date",
y = "Half-hour")
plotly::ggplotly(p)
Figure 4.2: Half hourly electricity import (current year, last 7 days)
plotDT[, dow := lubridate::wday(dv_date, label = TRUE)]
recent_dt <- plotDT[, .(sum_kWh_elec = sum(consumption_kWh)), keyby = .(dv_date, dow, dv_peakPeriod)]
daily_totals <- plotDT[, .(Total = sum(consumption_kWh)), keyby = .(dv_date, dow)]
t <- dcast(recent_dt, dv_date + dow ~ dv_peakPeriod, val.var = sum_kWh_elec)
## Using 'sum_kWh_elec' as value column. Use 'value.var' to override
# add totals
t <- t[daily_totals]
makeFlexTable(t, digits = 2,
cap = "Recent electricity use")
dv_date | dow | Early morning | Morning peak | Day time | Evening peak | Late evening | Total |
2022-09-13 | Tue | 2.69 | 0.75 | 0.92 | 2.85 | 1.66 | 8.87 |
2022-09-14 | Wed | 1.92 | 1.27 | 2.21 | 4.95 | 2.06 | 12.40 |
2022-09-15 | Thu | 2.00 | 0.72 | 1.25 | 3.76 | 1.60 | 9.32 |
2022-09-16 | Fri | 2.41 | 0.94 | 3.91 | 3.86 | 1.57 | 12.68 |
2022-09-17 | Sat | 1.86 | 0.89 | 1.75 | 1.84 | 2.89 | 9.22 |
2022-09-18 | Sun | 2.03 | 0.33 | 3.39 | 4.23 | 1.31 | 11.29 |
2022-09-19 | Mon | 1.56 | 0.96 | 2.47 | 3.50 | 1.38 | 9.87 |
plotDT <- elecCons_dt[, .(sum_kWh = sum(consumption_kWh),
mean_kWh = mean(consumption_kWh),
nObs = .N), keyby = .(dv_date, dv_weekend)]
ggplot2::ggplot(plotDT, aes(x = dv_date, y = mean_kWh,
colour = dv_weekend)) +
geom_point() +
geom_smooth() +
#facet_grid(dv_peakPeriod ~ .) +
scale_colour_viridis_d(name = "Weekend") +
labs(x = "Date",
y = "Mean kWh")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Figure 4.3: Daily electricity import (current year)
Figure 4.4 shows the mean daily kWh import with a smoothed curve for each period as defined below.
Early morning is effectively our baseload.
# check periods
t <- elecCons_dt[, .(min = min(dv_hms),
max = max(dv_hms)),
keyby = .(dv_peakPeriod)]
t
## dv_peakPeriod min max
## 1: Early morning 00:00:00 06:30:00
## 2: Morning peak 07:00:00 08:30:00
## 3: Day time 09:00:00 15:30:00
## 4: Evening peak 16:00:00 19:30:00
## 5: Late evening 20:00:00 23:30:00
plotDT <- elecCons_dt[, .(sum_kWh = sum(consumption_kWh),
mean_kWh = mean(consumption_kWh),
nObs = .N), keyby = .(dv_date, dv_peakPeriod, dv_weekend)]
ggplot2::ggplot(plotDT, aes(x = dv_date, y = mean_kWh,
colour = dv_peakPeriod)) +
geom_line() +
geom_smooth() +
#facet_grid(dv_peakPeriod ~ .) +
theme(legend.position = "bottom") +
guides(colour = guide_legend (ncol = 3)) +
scale_colour_viridis_d(name = "Peak period") +
labs(x = "Date",
y = "Mean kWh per period")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Figure 4.4: Daily electricity import by peak period (current year)
Not Octopus data - irregular reads by us. Load from .xslsx
Figure 4.5 compares our PV generation with overall grid import. These are mean values. We are occasionally exporting (Table 4.2 has zeros) but we do not (currently) have access to the export data.
library(openxlsx)
f <- "~/Dropbox/Home/houses/whiteHorseMews/2WhiteHorseMewsRunningCosts.xlsx"
pvGen <- openxlsx::read.xlsx(f,
sheet = "PV",
detectDates = TRUE)
pvGen_DT <- data.table::as.data.table(pvGen)
pvGen_DT[, meankWhGen := diff/n.days]
pvGen_DT[, dv_date := lubridate::as_datetime(Date)]
pvGen_DT[, month := lubridate::month(dv_date)]
monthlyPvGen_DT <- pvGen_DT[, .(hh_meankWh = mean(meankWhGen)/48 # input value is mean per day
),
keyby = .(month = lubridate::month(dv_date),
year = lubridate::year(dv_date))
]
monthlyPvGen_DT[, value := "PV generation"]
monthlyElec_DT <- elecCons_dt[, .(hh_meankWh = mean(consumption_kWh)),
keyby = .(month = lubridate::month(dv_date),
year = lubridate::year(dv_date))]
monthlyElec_DT[, value := "Grid import"]
plotDT <- rbind(monthlyPvGen_DT, monthlyElec_DT)
ggplot2::ggplot(plotDT, aes(x = month, y = hh_meankWh,
colour = value)) +
geom_line()+
scale_color_discrete(name = "Type")+
facet_grid(year ~ .)
## Warning: Removed 1 row(s) containing missing values (geom_path).
Figure 4.5: Compare grid import and PV gen
To do: model value of PV gen if we were to use all of it.
This will be a new MPAN but specified as export - although the url will still say consumption. We do not have this even though the PV is exporting on (some) days.
It may be that we only get this data if we sign up for the export tariff.
In theory our emissions from electricity use are zero because we are on a renewable-only tariff. But life is not so simple. We don’t have a private wire to a wind turbine so the electrons we import (stick with it) are as averagely green as all the rest.
We also con’t avoid the ‘Well To Tank’ emissions and those associated with transmission losses.
To further complicate things there are at least two different ways to estimate our emissions.
Does it matter? you cry. Well it might. If we’ve been able to ‘flex’ our usage in line with @theBakingForecast then who knows, maybe we’ll be concentrating our usage in times when the grid is actually drawing on more renewables.
So let’s take a look. We’ll do both the BEIS-based and NG-ESO based calculations to see. For now we’ll ignore the WTT and the T&D losses to keep the results comparable. We’ll come back to that later.
rmdParams$BEIS_elec_ci <- 0.21233
For the BEIS method, we’ll have to use the 2021 emissions factor as the 2022 value is not yet available.
For 2021 this is: 0.21233 Kg CO2e/kWh
For the NG-ESO method we use the NG-ESO half-hourly carbon intensity data that match to the half-hours in our electricity use dataset.
We save this locally, if it is older than 1 day we re-download.
# this needs to be more clever - only download dates we want from API?
esoF <- here::here("data", "latest_df_fuel_ckan.csv")
if(file.exists(paste0(esoF, ".gz"))){
message("We already have a version saved to: ", paste0(esoF, ".gz"))
message("Loading it...")
ngeso_dt_orig <- data.table::fread(paste0(esoF, ".gz"))
} else {
message("We don't already have a version, downloading and saving to: ", esoF)
ngeso_dt_orig <- data.table::fread("https://data.nationalgrideso.com/backend/dataset/88313ae5-94e4-4ddc-a790-593554d8c6b9/resource/f93d1835-75bc-43e5-84ad-12472b180a98/download/df_fuel_ckan.csv")
# nice dateTime
ngeso_dt_orig[, dv_start := lubridate::as_datetime(DATETIME)]
data.table::fwrite(ngeso_dt_orig, esoF)
dkUtils::gzipIt(esoF)
data.table::fwrite(ngeso_dt_orig, "~/Dropbox/data/UK_NGESO/genMix/latest_df_fuel_ckan.csv") # save locally for future re-use
dkUtils::gzipIt("~/Dropbox/data/UK_NGESO/genMix/latest_df_fuel_ckan.csv")
}
## We already have a version saved to: /Users/ben/repos/github/dataknut/octopusAPI/data/latest_df_fuel_ckan.csv.gz
## Loading it...
# if older than 1 day, reload
today <- lubridate::today()
lastNGESO <- as.Date(max(ngeso_dt_orig$dv_start))
if(today - lastNGESO > 1) {
# old data, reload
message("But the version we have dates from ", lastNGESO, " (",today - lastNGESO ," days ago), downloading latest...")
ngeso_dt_orig <- data.table::fread("https://data.nationalgrideso.com/backend/dataset/88313ae5-94e4-4ddc-a790-593554d8c6b9/resource/f93d1835-75bc-43e5-84ad-12472b180a98/download/df_fuel_ckan.csv")
# nice dateTime
ngeso_dt_orig[, dv_start := lubridate::as_datetime(DATETIME)]
data.table::fwrite(ngeso_dt_orig, esoF)
dkUtils::gzipIt(esoF)
}
# merge to usage data
setkey(ngeso_dt_orig, dv_start)
setkey(elecCons_dt, dv_start)
elecCons_dt <- ngeso_dt_orig[, .(dv_start, CARBON_INTENSITY, LOW_CARBON_perc, RENEWABLE_perc)][elecCons_dt] # keeps match to our electricity use
# there will be NAs if some datetimes are missing from ngeso_dt_orig
# we think renewable is wind + solar, low carbon includes nuclear
#ggplot2::ggplot(t, aes(x = RENEWABLE_perc, y = LOW_CARBON_perc)) +
# geom_point()
Mean half-hourly carbon intensity from the NG-ESO data for the data period was 0.1935 Kg CO2e/kWh. If this is substantially different to the BEIS 2021 value of 0.2123 Kg CO2e/kWh, we would expect emissions estimates using the NG-ESO factor to differ.
For context, Figure 4.6 summarises the mean half-hourly carbon intensity by month for the data period. We can clearly see that February 2022 was a very low carbon month… in fact it was a very windy month with 3 named storms.
elecCons_dt[, dv_month := lubridate::month(dv_date, label = TRUE)]
ggplot2::ggplot(elecCons_dt, aes(x = dv_month, y = CARBON_INTENSITY)) +
geom_violin(draw_quantiles = c(0.5)) +
#geom_jitter() +
#geom_boxplot() +
labs(x = "Month",
y = "Half-hourly carbon intensity",
caption = "Median drawn")
Figure 4.6: Monthly mean carbon intensity for the data period by month (NG-ESO data)
Figure 4.7 shows our half-hourly electricity kWh use vs halfhourly carbon intensity. Ideally we want a negative correlation showing that we use the most electricity when it is ‘greenest’ (carbon intensity is lowest). Doesn’t look too good, aye?
ggplot2::ggplot(elecCons_dt, aes(x = CARBON_INTENSITY, y = consumption_kWh, colour = RENEWABLE_perc)) +
geom_point() +
facet_wrap(. ~ dv_peakPeriod) +
geom_smooth() +
scale_color_continuous(name = "% renewables", low = "red", high = "green") +
theme(legend.position = "bottom") +
labs(y = "Halfhourly electricity kWh")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Figure 4.7: Own half-hourly electricity kWh vs NG-ESO halfhourly carbon intensity
What if we visualise using a box plot according to carbon intensity decile? So this means we divide the carbon intensity values into 10 equal groups - deciles. This is Figure 4.8. Doesn’t look too good either - our median usage (the horizontal bar in the boxes) seems to trend slightly upwards as we move to higher carbon intensity deciles.
# this will generate NAs if the CI data is missing for some of the (most recent) dateTimes
elecCons_dt[, CI_deciles := cut_number(CARBON_INTENSITY, n = 10)]
ggplot2::ggplot(elecCons_dt, aes(x = CI_deciles, y = consumption_kWh)) +
geom_boxplot() +
labs(y = "Halfhourly electricity kWh",
x = "Carbon intensity decile")
Figure 4.8: Boxing clever
So what if we just add up all our electricity kWh by carbon intensity decile? Do we use more low carbon kWh? This is Figure 4.9. Nah. The bakingforecast isn’t going to like us…
t <- elecCons_dt[, .(sumkWh = sum(consumption_kWh, na.rm = TRUE),
meankWh = mean(consumption_kWh, na.rm = TRUE)),
keyby = .(CI_deciles)]
ggplot2::ggplot(t, aes(x = CI_deciles, y = sumkWh)) +
geom_col() +
labs(y = "Sum kWh",
x = "Carbon intensity decile")
Figure 4.9: Sum of electricity kWh by carbon intensity decile
Out of interest, do our emissions values look very different if we apply the BEIS 2021 annual factor to our total electricity kWh to date compared to applying the NG-ESO half-hourly values?
elecCons_dt[, KgCO2_ngeso := consumption_kWh * (CARBON_INTENSITY/1000)] # convert to kg
t <- elecCons_dt[, .(sumkWh = sum(consumption_kWh),
sumKgCO2_ngeso = sum(KgCO2_ngeso, na.rm = TRUE))]
t[, sumKgCO2_beis := sumkWh * rmdParams$BEIS_elec_ci]
makeFlexTable(t, cap = "Comparing emissions estimation methods using electricity kWh to date")
sumkWh | sumKgCO2_ngeso | sumKgCO2_beis |
1,951 | 389 | 414 |
t <- elecCons_dt[, .(sumkWh = sum(consumption_kWh),
sumKgCO2_ngeso = sum(KgCO2_ngeso, na.rm = TRUE)),
keyby = .(dv_month)]
t[, sumKgCO2_beis := sumkWh * rmdParams$BEIS_elec_ci]
plotDT <- melt(t, id.vars = "dv_month")
ggplot2::ggplot(plotDT[ variable != "sumkWh",], aes(x = dv_month, y = value, fill = variable)) +
geom_col(position = "dodge") +
scale_color_discrete(name = "Method") +
labs(y = "Kg CO2",
x = "Month")
As we’d expect from the comparison of the values above, Table 4.4 suggests that it does. In fact our ‘in use’ NG-ESO based emissions are 65.28, 91.85, 87.3, 90.32, 88.2, 100.01, 106.32, 98.22 % of our BEIS-based emissions depending on the month in question.
If we compare the monthly values we can see the biggest difference was in February, a month we have already identified as being more ‘low carbon’ (see Figure 4.6).
We need to convert the gas consumption from m3 to kWh - see https://developer.octopus.energy/docs/api/#list-consumption-for-a-meter
gasM3TokWh <- 11.36
We use a multiplier of 11.36 kWh/m3 (https://www.theenergyshop.com/guides/how-to-convert-gas-units-to-kwh)
url <- paste0("https://api.octopus.energy/v1/gas-meter-points/",
apiParams$gas_mpan , "/",
"meters/",
apiParams$gas_serial, "/",
"consumption",
"?period_from=2022-01-01T00:00Z",
"&page_size=10000")
resp <- httr::GET(url = url, authenticate(user = apiParams$key, password = ""))
df <- jsonlite::parse_json(resp, simplifyVector = TRUE)
## No encoding supplied: defaulting to UTF-8.
gasCons_dt <- data.table::as.data.table(df$results)
gasCons_dt <- makeDerivedVars(gasCons_dt)
# gas 'consumption' is m3 - https://developer.octopus.energy/docs/api/#list-consumption-for-a-meter
# convert to kWh
gasCons_dt[, consumption_m3 := consumption]
gasCons_dt[, consumption_kWh := consumption * gasM3TokWh]
Note that this data starts later as we finally got the original un-registered smart meter replaced in February.
t <- gasCons_dt[, .(nDays = uniqueN(dv_date),
sumkWh = sum(consumption_kWh),
halfHourly_meankWh = mean(consumption_kWh)),
keyby = .(month = lubridate::month(dv_date))]
makeFlexTable(t, cap = "Monthly stats")
month | nDays | sumkWh | halfHourly_meankWh |
2 | 7 | 492 | 2 |
3 | 31 | 1,834 | 1 |
4 | 30 | 1,689 | 1 |
5 | 31 | 1,385 | 1 |
6 | 30 | 1,186 | 1 |
7 | 31 | 1,192 | 1 |
8 | 31 | 1,045 | 1 |
9 | 19 | 551 | 1 |
message("Total gas to date")
## Total gas to date
sum(t$sumkWh)
## [1] 9374.192
t <- gasCons_dt[, .(nDays = uniqueN(dv_date),
sumkWh = sum(consumption_kWh)),
keyby = .(dv_date)]
message("########")
## ########
gas_kWh <- 0.0719
message("Projected gas annual total @ ", gas_kWh,"p")
## Projected gas annual total @ 0.0719p
annualTotal <- mean(t$sumkWh)*365
annualTotal
## [1] 16293.24
message("Projected annual gas cost @ ", gas_kWh,"p - including standing charge")
## Projected annual gas cost @ 0.0719p - including standing charge
annualCost <- (annualTotal*gas_kWh)+(365*0.261) # standing charge
annualCost
## [1] 1266.749
message("Projected mean monthly gas cost @ ", gas_kWh,"p - including standing charge")
## Projected mean monthly gas cost @ 0.0719p - including standing charge
annualCost/12
## [1] 105.5624
message("########")
## ########
gas_kWh <- 0.103
message("Projected gas annual total @ ", gas_kWh,"p")
## Projected gas annual total @ 0.103p
annualTotal <- mean(t$sumkWh)*365
annualTotal
## [1] 16293.24
message("Projected annual gas cost @ ", gas_kWh,"p - including standing charge")
## Projected annual gas cost @ 0.103p - including standing charge
annualCost <- (annualTotal*gas_kWh)+(365*0.261) # standing charge
annualCost
## [1] 1773.469
message("Projected mean monthly gas cost @ ", gas_kWh,"p - including standing charge")
## Projected mean monthly gas cost @ 0.103p - including standing charge
annualCost/12
## [1] 147.7891
Figure 4.10 shows half-hourly gas import (‘consumption’) for the current year. The power cuts are even easier to see here. Interestingly the pattern after the gas boiler was serviced in June is more varied. What did he change?
ggplot2::ggplot(gasCons_dt, aes(x = dv_date, y = dv_hms, fill = consumption_kWh)) +
geom_tile() +
theme(legend.position = "bottom") +
scale_fill_viridis_c(name = "Gas import (kWh)") +
labs(x = "Date",
y = "Half-hour")
Figure 4.10: Half-hourly gas consumption (current year)
Repeat but with just the last 7 days of data - useful for checking recent appliance use and offspring effects.
today <- lubridate::today()
plotDT <- gasCons_dt[dv_date >= today - 7]
p <- ggplot2::ggplot(plotDT, aes(x = dv_date, y = dv_hms, fill = consumption_kWh)) +
geom_tile() +
theme(legend.position = "bottom") +
scale_fill_viridis_c(name = "Gas import (kWh)") +
labs(x = "Date",
y = "Half-hour")
plotly::ggplotly(p)
Figure 4.11: Half hourly gas import (current year, last 7 days)
plotDT[, dow := lubridate::wday(dv_date, label = TRUE)]
recent_dt <- plotDT[, .(sum_kWh_elec = sum(consumption_kWh)), keyby = .(dv_date, dow, dv_peakPeriod)]
daily_totals <- plotDT[, .(Total = sum(consumption_kWh)), keyby = .(dv_date, dow)]
t <- dcast(recent_dt, dv_date + dow ~ dv_peakPeriod, val.var = sum_kWh_elec)
## Using 'sum_kWh_elec' as value column. Use 'value.var' to override
# add totals
t <- t[daily_totals]
makeFlexTable(t, digits = 2,
cap = "Recent gas use")
dv_date | dow | Early morning | Morning peak | Day time | Evening peak | Late evening | Total |
2022-09-13 | Tue | 12.28 | 0.00 | 0.00 | 13.84 | 0.00 | 26.12 |
2022-09-14 | Wed | 12.25 | 0.00 | 0.00 | 16.73 | 0.00 | 28.98 |
2022-09-15 | Thu | 17.37 | 0.25 | 0.00 | 16.49 | 0.00 | 34.11 |
2022-09-16 | Fri | 13.35 | 0.00 | 0.00 | 19.28 | 0.00 | 32.63 |
2022-09-17 | Sat | 13.36 | 0.00 | 0.00 | 18.83 | 0.00 | 32.19 |
2022-09-18 | Sun | 18.56 | 0.17 | 0.00 | 16.18 | 0.00 | 34.91 |
2022-09-19 | Mon | 13.72 | 0.00 | 0.00 | 16.28 | 0.00 | 30.00 |
Figure 4.12 shows the mean daily kWh import with a smoothed curve.
To do: mark weekends etc
plotDT <- gasCons_dt[, .(sum_kWh = sum(consumption_kWh),
mean_kWh = mean(consumption_kWh),
nObs = .N), keyby = .(dv_date, dv_weekend)]
ggplot2::ggplot(plotDT, aes(x = dv_date, y = sum_kWh,
colour = dv_weekend)) +
geom_point() +
geom_smooth() +
theme(legend.position = "bottom") +
guides(colour = guide_legend (ncol = 3)) +
scale_colour_viridis_d(name = "Weekend") +
labs(x = "Date",
y = "Sum kWh per day")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Figure 4.12: Daily gas consumption (current year)
Repeat for mean
ggplot2::ggplot(plotDT, aes(x = dv_date, y = mean_kWh,
colour = dv_weekend)) +
geom_point() +
geom_smooth() +
theme(legend.position = "bottom") +
guides(colour = guide_legend (ncol = 3)) +
scale_colour_viridis_d(name = "Weekend") +
labs(x = "Date",
y = "Mean kWh per day")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Figure 4.13 shows the mean daily kWh import with a smoothed curve by period of the day.
To do: mark weekends etc
plotDT <- gasCons_dt[, .(sum_kWh = sum(consumption_kWh),
mean_kWh = mean(consumption_kWh),
nObs = .N), keyby = .(dv_date, dv_peakPeriod, dv_weekend)]
ggplot2::ggplot(plotDT, aes(x = dv_date, y = mean_kWh,
colour = dv_peakPeriod)) +
geom_line() +
geom_smooth() +
#facet_grid(dv_peakPeriod ~ .) +
theme(legend.position = "bottom") +
guides(colour = guide_legend (ncol = 3)) +
scale_colour_viridis_d(name = "Peak period") +
labs(x = "Date",
y = "Mean kWh per period")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Figure 4.13: Daily gas consumption by peak period (current year)
This is much more simple. We can only apply the BEIS 2021 value as there are no time-varying emissions factors for gas.
rmdParams$BEIS_gas_ci <- 0.20297
As before, for the BEIS method we’ll have to use the 2021 emissions factor as the 2022 value is not yet available.
For 2021 this is: 0.20297 Kg CO2e/kWh
gasCons_dt[, KgCO2_beis := consumption_kWh * rmdParams$BEIS_gas_ci]
t <- gasCons_dt[, .(sumkWh = sum(consumption_kWh),
sumKgCO2_beis = sum(KgCO2_beis))]
makeFlexTable(t, cap = "Emissions estimation using gas kWh to date")
sumkWh | sumKgCO2_beis |
9,374 | 1,903 |
Use skmir::skim()to summarise.
skimr::skim(elecCons_dt)
| Name | elecCons_dt |
| Number of rows | 10000 |
| Number of columns | 15 |
| Key | dv_start |
| _______________________ | |
| Column type frequency: | |
| character | 3 |
| Date | 1 |
| difftime | 1 |
| factor | 3 |
| numeric | 6 |
| POSIXct | 1 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| interval_start | 0 | 1 | 20 | 25 | 0 | 10000 | 0 |
| interval_end | 0 | 1 | 20 | 25 | 0 | 10000 | 0 |
| dv_weekend | 0 | 1 | 6 | 8 | 0 | 3 | 0 |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dv_date | 0 | 1 | 2022-02-23 | 2022-09-19 | 2022-06-07 | 209 |
Variable type: difftime
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dv_hms | 0 | 1 | 0 secs | 84600 secs | 43200 secs | 48 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| dv_peakPeriod | 0 | 1 | FALSE | 5 | Ear: 2912, Day: 2912, Eve: 1672, Lat: 1672 |
| dv_month | 0 | 1 | TRUE | 8 | Mar: 1488, May: 1488, Jul: 1488, Aug: 1488 |
| CI_deciles | 0 | 1 | FALSE | 10 | (10: 1024, [39: 1014, (24: 1001, (23: 1000 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| CARBON_INTENSITY | 0 | 1 | 193.52 | 58.71 | 39.0 | 153.00 | 205.00 | 239.00 | 323.00 | ▂▃▆▇▂ |
| LOW_CARBON_perc | 0 | 1 | 50.94 | 13.71 | 22.0 | 40.00 | 48.80 | 60.30 | 87.70 | ▂▇▆▃▂ |
| RENEWABLE_perc | 0 | 1 | 28.79 | 14.64 | 2.2 | 17.00 | 27.30 | 39.00 | 70.20 | ▆▇▆▃▁ |
| consumption | 0 | 1 | 0.20 | 0.20 | 0.0 | 0.08 | 0.14 | 0.23 | 1.79 | ▇▁▁▁▁ |
| consumption_kWh | 0 | 1 | 0.20 | 0.20 | 0.0 | 0.08 | 0.14 | 0.23 | 1.79 | ▇▁▁▁▁ |
| KgCO2_ngeso | 0 | 1 | 0.04 | 0.04 | 0.0 | 0.01 | 0.03 | 0.05 | 0.42 | ▇▁▁▁▁ |
Variable type: POSIXct
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dv_start | 0 | 1 | 2022-02-23 16:00:00 | 2022-09-19 23:30:00 | 2022-06-07 19:45:00 | 10000 |
Use skmir::skim()to summarise.
skimr::skim(gasCons_dt)
| Name | gasCons_dt |
| Number of rows | 10000 |
| Number of columns | 11 |
| Key | NULL |
| _______________________ | |
| Column type frequency: | |
| character | 3 |
| Date | 1 |
| difftime | 1 |
| factor | 1 |
| numeric | 4 |
| POSIXct | 1 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| interval_start | 0 | 1 | 20 | 25 | 0 | 10000 | 0 |
| interval_end | 0 | 1 | 20 | 25 | 0 | 10000 | 0 |
| dv_weekend | 0 | 1 | 6 | 8 | 0 | 3 | 0 |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dv_date | 0 | 1 | 2022-02-22 | 2022-09-19 | 2022-06-07 | 210 |
Variable type: difftime
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dv_hms | 0 | 1 | 0 secs | 84600 secs | 43200 secs | 48 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| dv_peakPeriod | 0 | 1 | FALSE | 5 | Day: 2913, Ear: 2912, Eve: 1672, Lat: 1671 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| consumption | 0 | 1 | 0.08 | 0.14 | 0 | 0 | 0 | 0.12 | 1.19 | ▇▁▁▁▁ |
| consumption_m3 | 0 | 1 | 0.08 | 0.14 | 0 | 0 | 0 | 0.12 | 1.19 | ▇▁▁▁▁ |
| consumption_kWh | 0 | 1 | 0.94 | 1.65 | 0 | 0 | 0 | 1.31 | 13.53 | ▇▁▁▁▁ |
| KgCO2_beis | 0 | 1 | 0.19 | 0.33 | 0 | 0 | 0 | 0.27 | 2.75 | ▇▁▁▁▁ |
Variable type: POSIXct
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dv_start | 0 | 1 | 2022-02-22 15:30:00 | 2022-09-19 23:00:00 | 2022-06-07 19:15:00 | 10000 |
Figure 5.1 shows the NG-ESO half-hourly carbon intensity over time for the data period as context.
ggplot2::ggplot(elecCons_dt, aes(x = dv_date, y = dv_hms, fill = CARBON_INTENSITY)) +
geom_tile() +
scale_fill_continuous(name = "Carbon intensity", low = "green", high = "red") +
labs(x = "Date",
y = "Time of day",
caption = "Source: NG-ESO (https://data.nationalgrideso.com/carbon-intensity1/historic-generation-mix)")
Figure 5.1: Half-hourly carbon intensity over time for the data period